function [thetastar,fval,exitflag,output] = dreumScript(Ncpu)

% Filename: dreumScript.m
% Author: Carson Reeling, Purdue University and Valentin Verdier, University of North Carolina - Chapel Hill
% Date last modified: Jan. 2020
% Description: Outer script for running MLE routine accompanying Reeling et
%   al. "Valuing Goods Accessed via Dynamic Lottery." Journal of the
%   Association of Environmental and Resource Economists.

%*********************************************
%To run on Dogwood at UNC (replace with own cluster initialization code):
%*********************************************

%Delete files in /21dayscratch/scr/v/v/vverdier/R2018a
%Submit job with "sbatch run.sl"

c = parcluster('dogwood R2019a');
delete(c.Jobs)
% 
c.AdditionalProperties.QueueName = 'skylake';
c.AdditionalProperties.WallTime = '6-00:00:00';
c.AdditionalProperties.AdditionalSubmitArgs = '-o out.parjob.%J -e err.parjob.%J';
% %--reservation=vverdier_43';
% %--qos workshop_access
% %--mem-per-cpu=14G;
% 
p = c.parpool(Ncpu);
% 

%*********************************************
%END
%*********************************************

global muHistory gradHistory LLHistory P23maxHistory P0maxHistory % Initialize variable to store intermediate values of theta during optimization 
gradHistory = [];
LLHistory = [];
P23maxHistory = [];
P0maxHistory = [];

% Set parameters
K = 10;                         % K - 1 = max number of preference points
H = 23;                         % Number of hunts, including no-apply option (except pp only option)
tau = 3;                        % Number of different "types" of hunters

rng(1234)

%Initial parameter vector
theta0 = [-.2, 3, 2, .03, .8, .2*[5, 3, 3, 7, 5, 4, 7, 4, 4, 6, 4, 3, 7, 4, 3, 9, 6, 5, 5, 9, 9, 5], .22, 5];
muHistory = theta0;

tolInner = 1e-8;               	% Tolerance for inner loop
rho=0.975;						% Discount factor

[choiceFull,Nu,pp,ppFull,tc,~,tcidx,~] = dreumDataProc3(1);

[chunksize,tcSplit] = dreumDataSplit2(Ncpu,Nu,tc);
         
% Load and prepare data
% choice................(N x 2) vector of applicants' application choices 
% data..................(N x dim) matrix of applicant data 
% Nu....................Number of unique applicants (defined by number of
%                       unique travel cost combinations)
% pp....................(N x 2) vector of applicants' preference point stock
% tc....................contains all covariates of the model, i.e. MC, OC, and applicant characteristics (income, age, gender, geographical region 
%						indicator variables, nearby distance indicator)
% tcidx.................(N x 1) vector of indices relating applicants'
%                       travel cost combinations to row in tc
 
[phi,PMarkov] = dreumProb2(H,K);  

% 	Set up division of applicants for parallelization (each part of parfor handles a fraction of all applicants)


% Calculate success probabilities and Markov transition matrices-----------

% phi...................(H+1 x K x 1) array of success probabilities for
%                       1st round for each preference point level
% PMarkov...............(K x K x H*H+1) array of Markov transition matrices
%                       relating the probability of transititioning from p 
%                       to p+1 preference points, conditional on hunt choice
        

% Estimate theta via maximum likelihood------------------------------------
% Define objective function
obj = @(x)dreumLL2_wg(choiceFull,ppFull,tcSplit,H,K,phi,PMarkov,rho,tau,...
    x,tolInner,Ncpu,Nu,tcidx,chunksize); % Specify objective function

            
% Set optimization options
% options = optimset('Display','iter','GradObj','on','OutputFcn',@outfun);
options = optimoptions(@fmincon,'Display','iter','SpecifyObjectiveGradient',true,'OutputFcn',@outfun);


% Set constraints
% Key: [mu      eta1 eta2 pi1 pi2 chi1-22         alpha s (scale parameter on exponential distribution)] 
ub =   [0       30  30    1   1   20*ones(1,22)  .3    Inf];    % Upper bound on parameter values
lb =   [-Inf    0    0      0   0   -20*ones(1,22) 0     0];      % Lower bound on parameter values
A =    [0       0    0    1   1   zeros(1,22)    0     0; 
        0       -1   1    0   0   zeros(1,22)    0     0;
        0       -1    0    0   0   zeros(1,22)    0     0];
b = [1 0 0]';


% Call minimization routine
% [thetastar,negLLstar,exitflag,output] = knitromatlab(obj,theta0,A,b,[],[],lb,ub,[],[],options)
[thetastar,negLLstar,exitflag,output] = fmincon(obj,theta0,A,b,[],[],lb,ub,[],options);

save('output.mat'); % Save estimation results

% delete(gcp('nocreate'));